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Abstract. In this paper we propose two different primal-dual splitting algorithms 
for solving inclusions involving mixtures of composite and parallel-sum type monotone 
operators which rely on an inexact Douglas-Rachford splitting method, however applied 
in different underlying Hilbert spaces. Most importantly, the algorithms allow to process 
the bounded linear operators and the set-valued operators occurring in the formulation 
of the monotone inclusion problem separately at each iteration, the latter being individ- 
ually accessed via their resolvents. The performances of the primal-dual algorithms are 
emphasized via some numerical experiments on location and image deblurring problems. 
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1 Introduction and preliminaries 

In applied mathematics, a wide range of convex optimization problems such as single- or 
multifacility location problems, support vector machine problems for classification and 
regression, portfolio optimization problems as well as signal and image processing prob- 
lems, all of them likely possessing nondifferentiable convex objectives, can be reduced 
to the solving of inclusions involving mixtures of monotone set- valued operators. 

In this article we propose two different primal-dual iterative error-tolerant methods 
for solving inclusions with mixtures of composite and parallel-sum type monotone op- 
erators. Both algorithms rely on the inexact Douglas-Rachford algorithm (cf. [8,9]), 
but still differ clearly from each other. An important feature of the two approaches 
and, simultaneously, an advantage over many existing methods is their capability of 
processing the set- valued operators separately via their resolvents, while the bounded 
linear operators are accessed via explicit forward steps on their own or on their adjoints. 
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The resolvents of the maximally monotone operators are not always available in closed 
form expressions, fact which motivates the inexact versions of the algorithms, where 
implementation errors in the shape of summable sequences are allowed. 

The methods in this article are also perfectly parallelizable since the majority of their 
steps can be executed independently. Furthermore, when applied to subdifferential 
operators of proper, convex and lower semicontinuous functions, the solving of the 
monotone inclusion problems is, under appropriate qualification conditions (cf. [2,4]), 
equivalent with finding optimal solutions to primal-dual pairs of convex optimization 
problems. The considered formulation also captures various different types of primal 
convex optimization problems and corresponding conjugate duals appearing in wide 
ranges of applications. The resolvents of subdifferentials of proper, convex and lower 
semicontinuous functions are the proximal point mappings of these and are known to 
take in a lot of situations shape of closed form expressions. 

Recent research (see [3,5,6,10,14]) has shown that structured problems dealing with 
monotone inclusions can be efficiently solved via primal-dual splitting approaches. In 
[10], the problem involving sums of set-valued, composed, Lipschitzian and parallel-sum 
type monotone operators was decomposed and solved via an inexact Tseng algorithm 
having foward-backward-forward characteristics in a product Hilbert space. On the 
other hand, in [14], instead of Lipschitzian operators, the author has assumed cocoercive 
operators and solved the resulting problem with an inexact forward- backward algorithm. 
Thus, our methods can be seen as a continuation of these ideas, this time by making 
use of the inexact Douglas-Rachford method. Another primal-dual method relying on 
the same fundamental splitting algorithm is considered in [11] in the context of solving 
minimization problems having as objective the sum of two proper, convex and lower 
semicontinuous functions, one of them being composed with a bounded linear operator. 

Due to the nature of Douglas-Rachford splitting, we will neither assume Lipschitz 
continuity nor cocoercivity for any of the operators present in the formulation of the 
monotone inclusion problem. The resulting drawback of not having operators which 
can be processed explicitly via forward steps is compensated by the advantage of allow- 
ing general maximal monotone operators in the parallel-sums, fact which relaxes the 
working hypotheses in [10, 14]. 

The article is organized as follows. In the remaining of this section we introduce the 
framework we work within and some necessary notations. The splitting algorithms and 
corresponding weak and strong convergence statements are subject of Section 2 while 
Section 3 is concerned with the application of the two methods to convex minimization 
problems. Finally, in Section 4 we make some numerical experiments and evaluate the 
obtained results. 

We are considering the real Hilbert spaces % and Qi endowed with the inner prod- 
uct {■,-)y^ and and associated norm ||-||^ = (•, •)^ and = = 
l,...,m, respectively. The symbols ^ and — ^ denote weak and strong convergence, 
respectively, M-|__|- denotes the set of strictly positive real numbers and IR+ = M-|_-|- U {0}. 
By -6(0, r) we denote the closed ball with center and radius r G For a function 
M = M U {±00} we denote by dom/ := {x £ Ti : f{x) < +00} its effective 
domain and call / proper if dom / / and f{x) > —00 for all x £ Ti. Let be 

r('H) := {/ : — ^ M : / is proper, convex and lower semicontinuous}. 
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The conjugate function of / is /* : — )■ M, f*{p) = sup{{p,x) — f{x) : x £ H} for 
all p £ Ti and, if / G ^{'H), then /* G r('H), as well. The (convex) subdifferential of 
/ : -H ^ M at X G is the set df{x) = {p £ U : f{y) - f{x) > {p,y- x) Vy G n}, if 
f{x) G M, and is taken to be the empty set, otherwise. For a linear continuous operator 
Li : H ^ Gi, the operator L* : Gi ^ H, defined via {LiX,y) = {x,L*y) for all x G 
and all y £ Gi, denotes its adjoint, for i G {1, . . . , m}. 

Having two functions /, 5 : H — > M, their infimal convolution is defined by fOg : 

— M, {fDg){x) = infj^g-^ {f{y) + — y)} for all x £ H, being a convex function 
when / and g are convex. 

Let M : "H — )■ 2^ be a set-valued operator. We denote by zer M = {x £ Ti : £ Mx} 
its set of zeros, hj &xM = {x £ 7i : x £ Mx} its set of fixed points, by gra M = {{x,u) £ 
7i X Ti : u £ Mx} its graph and by ranM = {u £ 7i : 3x £ H, u £ Mx} its range. 
The inverse of M is : U ^ 2^ , u <^ {x £ U : u £ Mx}. We say that the 

operator M is monotone ii {x — y,u — v) > for all (y, w) G graM and it is 

said to be maximally monotone if there exists no monotone operator M' : % ^ 2^ 
such that graM' properly contains graM. The operator M is said to be uniformly 
monotone with modulus (/>m : 1^+ [0,-|-oo] if (j)M is increasing, vanishes only at 0, 
and {x — y,u — v) > 4>m {\\x — y\\) for all (x, u), {y, v) £ gra M. 

The resolvent and the reflected resolvent of an operator M : H — >■ 2^ are 

JAf = (Id + M)-^ and Rm = 2Jm - Id, 

respectively, the operator Id denoting the identity on the underlying Hilbert space. 
When M is maximally monotone, its resolvent (and, consequently, its reflected resol- 
vent) is a single-valued operator and, by [1, Proposition 23.18], we have for 7 G 

Id = + 7-^7-1 A/-1 ° 7"4d. (1.1) 

Moreover, for / G ^{"H) and 7 G the subdifferential ^(7/) is maximally monotone 
(cf. [15, Theorem 3.2.8]) and it holds J-yQf = (Id -|- 79/)^"'^ = Prox^j. Here, Prox^j(x) 
denotes the proximal point of 7/ at x G representing the unique optimal solution of 
the optimization problem 

mf |7/(y) + ^||y-x||2|. (1.2) 

In this particular situation (1.1) becomes Moreau's decomposition formula 

Id = Prox^j -I-7 Prox^-i o7^^Id. (1-3) 

When $7 C is a nonempty, convex and closed set, the function 6n : 7i ^ M, defined 
by Sq{x) = for X G n and 5q(x) = -|-oo, otherwise, denotes the indicator function of 
the set J7. For each 7 > the proximal point of at x G H is nothing else than 

Prox^5j^(x) = Prox5j,(x) = Vnix) = argmin - x||^, 

yen ^ 

which is the projection of x on 
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The sum and the parallel sum of two set- valued operators Mi, M2 : 7^ — >■ 2^ are 
defined as Mi + M2 : H ^ 2'^, (Mi + M2)(x) = Mi(x) + M2{x) Vrc e "H and 

Ml □ M2 : -H ^ 2^, Ml □ M2 = (Mf ^ + Ma"^) , 

respectively. If Mi and M2 are monotone, than Mi + M2 and Mi □ M2 are monotone, 
too, however, if Mi and M2 are maximally monotone, this property is in general not 
true neither for Mi + M2 nor for Mi □ M2 (see [2] ) . 



2 Algorithms and convergence results 

Within this section we provide two algorithms together with weak and strong conver- 
gence results for the following primal-dual pair of monotone inclusion problems. 

Problem 2.1. Let A : T-L ^ 2^ he a maximally monotone operator and z ^ T-L. 
Furthermore, for every i G {1, . . . ,m}, let ri G Qi, Bi : Qi ^ 2^' and Di : Qi ^ 2^' be 
maximally monotone operators and Li : % ^ Qi a, nonzero bounded linear operator. 
The problem is to solve the primal inclusion 

m 

find X G -H such that z & Ax + L*{Bi □ Di){LiX - n) (2.1) 

i=l 

together with the dual inclusion 



find vi £ Qi, . . . ,Vm ^ Gm such that {3x £ H) 



_z-ET=iL*v,eAx 

Vi e{BiUDi){LiX - n), i = l,...,m. 

(2.2) 



We say that (x, f 1, . . . , Vm) & Ti x Qi . . . x Qm is a primal-dual solution to Problem 
2.1, if 

m 

z — ^ L*Vi G Ax and Vi G {Bi □ Di){LiX — r^), i = 1, . . . , m. (2-3) 

i=l 

If (x, ui, . . . ,Vrn) £ Ti X Qi . . . X Q^n is a primal-dual solution to Problem 2.1, then x is 
a solution to (2.1) and (vi, . . . ,Vm) is a solution to (2.2). Notice also that 

m 

X solves (2.1) 4^ z - ^ L* {BiD Di){LiX - r^) e Ax 4^ 
1=1 

^ z-YT=i L*Vi G Ax, 
Vi e {BiUDi){LiX - Ti), i = l,...,m. 



3vi e Qi, . . . ,Vm & Gm such that 



Thus, if X is a solution to (2.1), then there exists (tJi, . . . ,Vm) G ^1 x . . . Gm such that 
{x,vi, . . . , Vm) is a primal-dual solution to Problem 2.1 and if (iJi, . . . , Vm) € Gix . . . Gm 
is a solution to (2.2), then there exists x G "H such that (x, ui, . . . ,Vm) is a primal-dual 
solution to Problem 2.1. 
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Example 2.1. In Problem 2.1, set m = 1, z = and ri = 0, let Qi = Q, Bi = B, 
Li = L, Di : g ^ 2^, Di{0) = G and Di{v) = \/v e G \ {0}, and ^ : H ^ 2^ and 
B : g ^ 2^ he the convex subdifferentials of the functions / G r('H) and g G r(^), 
respectively. Then, under appropriate qualification conditions (see [2,4]), to solve the 
primal inclusion problem (2.1) is equivalent to solve the optimization problem 

inf {/(x)+5(Lx)}, 

while to solve the dual inclusion problem (2.2) is nothing else than to solve its Fenchel- 
dual problem 

sup{-n-L*v)-g*iv)}. 

For more primal-dual pairs of convex optimization problems which are particular 
instances of (2.1)-(2.2) we refer to [10,14]. 



2.1 A first primal-dual algorithm 

The first iterative scheme we propose in this paper and which we describe as follows 
has the particularity that it accesses the resolvents of A, B^^ and D~^,i = 1, m, and 
processes each operator Li and its adjoint L*,i = l,...,m two times. 

Algorithm 2.1. 

Let xq £ Ti, (vifi, . . . , Vmfl) £ GiX . . .x Qm and r and dj, i = 1, m, be strictly positive 
real numbers such that 

m 

T^a.WLif < 4. 
1=1 

Furthermore, let (An)n>o be a sequence in (0, 2), (a„)n>o a sequence in 7i, {bi^n)n>o and 
{di,n)n>o sequences in Gi for all i = 1, . . . , m and set 

Pl,n = JtA {Xn - \ Ya=\ ^'>^i,n + Tz) + a„ 

For i = 1, . . . ,m 

* i 

W2,i,n = 2p2,j,n " Vi^n 
Xn+1 = Xn -\- ^n{zi^n ~ Pl,n) 

For i = 1 , . . . , m 

Z2,i,n = J„ r)-^ {W2,i,n + ^Li{2zi^n - t«l,n)) + ^i.n 
Vi,n+1 = Vi,n + An(^2,i,n — P2,i,n)- 



(Vn > 0) 



(2.4) 



Theorem 2.1. For Problem 2.1 assume that 

Cm 
A + Y,L*{B,UDi){Li 
i=i 

and consider the sequences generated by Algorithm 2.1. 



(2.5) 
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(i) If 

+00 +00 

X! < +°0> X! '^n{\\di^n\\g, + \\h,n\\Q^) < +00,i = 1, . . . ,m, 
n=0 n=0 

and I]n:S A„(2 - An) = +00, then 

(a) [xni vi^n, ■ ■ ■ 1 Vm,n)n>o Converges weakly to a point {x,vi, . . . ,Vm) £ x Qi x 
... X Qm such that, when setting 

Pi = -^tA - ^ ^ L*Vi + rz j , 

and P2^i = Ja,B-^ {^i + -x) - atr,^ , i = l, ■■■,m, 

the element {pi,P2 i; • • • )P2 m) ^'-^ 0, primal-dual solution to Problem 2.1. 

(b) A„(zi,„ - pi,n) (n +00) and Xn{z2,i,n - P2,i,n) (n +00) for 
i = 1, m. 

(c) whenever % and Qi, i = l,...,m, are finite- dimensional Hilbert spaces, 
0"n — >■ (n — )■ +00) and bi^n — )• (n — )■ +00) for i = then 
{pi,mP2,i,n, ■ ■ ■ ,P2,m,n)n>o convcrges strongly to a primal-dual solution of 
Problem 2.1. 

(ii) If 

+00 +00 

f lid; r.. lie- + II&7; -n lL, ) < +00, i = 1 ,m, 

n>0 



^\\an\\H< +00, ^{\\di^n\\g, + \\bi,n\\g,) < +oo,i = I, . . . ,m, infA„>0 

n=0 n=0 



and A and B- = 1, ...,m, are uniformly monotone, 

then {pi,n,P2,i,n, ■ ■ ■ ,P2,m,n)n>o convcrgcs strongly to the unique primal-dual solu- 
tion of Problem 2.1. 

Proof. Consider the Hilbert space Q = Gi x . . . x Qm endowed with inner product and 
associated norm defined, for v = {vi, . . . , Vm), Q = {li, ■ ■ ■ 1 Qm) G Qi as 



{'"^(l)g = ^{vi,Qi)g^ and ||i>||g 



i=l 



.EMh (2-6) 
\ i=i 



respectively. Furthermore, consider the Hilbert space K, = % x Q endowed with inner 
product and associated norm defined, for {x,v),{y,q) S /C, as 



{{x,v),{y,q)),^ = {x,y)y^ + {v,q)g and \\{x,v)\\}c = ylkH^ + ||-w||^, (2.7) 
respectively. Consider the set-valued operator 

,vi, . . . , Vm) ^ {—z + Ax, ri + Vi, . . ■ , r^, + B 
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which is maximally monotone, since A and Bi,i = l,...,m, are maximally monotone 
(cf. [1, Proposition 20.22 and Proposition 20.23]) and the bounded linear operator 

S :K^K, {x,Vi,..., Vm) ^ --^1^' • • • ' -LmX^ , 

which proves to be skew (i. e. S* = —S) and hence maximally monotone (cf. [1, Example 
20.30]). Further, consider the set- valued operator 

Q:K^2^, {x,vi,...,Vm)^ {{),D^\i, . . . ,D^\. 



which is maximally monotone, as well, since Di is maximally monotone for i = 1, m. 
Therefore, since domS = /C, both + Q and + M are maximally monotone 
(cf. [1, Corollary 24.4(i)]). On the other hand, according to [10, Eq. (3.12)], it holds 
(2.5) <^ zer (M + 5 + Q) / 0, while [10, Eq. (3.21) and (3.22)] yield 

(x, f 1, . . . , Vm) £ zer [M + S + Q) =^(a;, ui, . . . , Vm) is a primal-dual 

solution to Problem 2.1. (2.8) 

Finally, we introduce the bounded linear operator 

V-.K^K, {x,Vi,...,Vm) ^ {- -\^L*Vi, — -]-Lix,. . . , — -]-Lmx\ . 

y ^ ~[ cri 2 CTm 2 J 

It is a simple calculation to prove that V is self-adjoint, i.e. V* = V. Furthermore, 
the operator V is /j-strongly positive for 



' 2 



1 nr 112 ■ / 



1 1 1 



T Gi On 



which is a positive real number due to the assumption 



m 

T^ai||L,f<4 (2. 

2=1 



made in Algorithm 2.1. Indeed, using that 2a6 < aa} -|- ^ for any a, 6 G M and any 
a E it yields for each i = 1, . . . ,m 

ollr nil II II II ^ II ||2 I V'^ 'l2iLl '^ill-^'iP n ||2 i n\ 

2||^i|||k||wlbi|l6. < / llr ^ll^ll^ + ll^ilb, (2.10) 
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and, consequently, for each x = {x,vi, . . . , Vm) S /C, it follows that 



> 



+ E 

i=l 
m 

+ E 

i=l 



\Vi\ 



1=1 

m 

E 



x\\n\\vi\\gi 




> 



^\ ~{ I It ai am 



\x\ 



(2.11) 



Since V is /9-strongly positive, we have cl(ran V) = ran V (cf. [1, Fact 2.19]), zeiV = 
{0} and, as (ran V)-*- = zer V* = zeiV = {0} (see, for instance, [1, Fact 2.18]), it holds 
ran V = JC. Consequently, exists and < ^. 

The algorithmic scheme (2.4) is equivalent to 



(Vn > 0) 



r 2 ^i=l ^i^i,n & A{Pl,n — ^r. 



Z 



For i = 1, . . . ,m 

- ^LiiXn - Pl,n) G -^LiPi^n + B-\p2,i,n " 6,,„) + 

W2,i,n = 2p2,i,n " Wj,n 

^— 2 ^i=l ^i^'i,i,n = U 

2;n+l = a^n + A„(zi,„ - Pl,n) 

For i = 1, . . . ,m 

"'^'''"g-/^'''" ~ 5-^j(^^l,n — Zi^n) S — 5-^i2;i,n + '^{z2,i,n — di,n) 
Vi,n+1 = Vi^n + Ki{z2,i,n — P2,i,n)- 



(2.12) 



We introduce for every n > the following notations: 



and < 



•^n — (a^n) ^l,n) • • • 5 ^m,n) 

= (Pl,n,P2,l,n, • • ■,P2,m,n) 

— (^l,ra; ^2,l,ri) • • • ) -2'2,m,?i) 



The scheme (2.12) can equivalently be written in the form 



dn = 


(0,(ii,„, . . 


) dm,n) 


dC^ — 


(0,^,.. 


' fr™ ) 


bn = 


(o-n) bi^iii • 


• ) bui fi 


[ K = 


(an h,n 
V T ' O-i ' • 





(Vn > 0) 



V{x^ - yj G (^5 + m) (y„ - bn) + \SK - hi 

Wn = 2y„ - Xn 

V{Wn - Zn) G [\S + Q) {Zn - dn) + \Sdn " < 
— "^n ~l~ i,Zn Vn) • 



(2.13) 



(2.14) 



8 



We set for every n > 



e°- = V 



-1 



S + v]dn-dl 



(2.15) 



Next we introduce the Hilbert space /Cy with inner product and norm respectively 
defined, for a;, y G /C, as 



y)Kv = ^y)K and \\x\\ic^ = ^ {x, Vx)^, (2.16) 
respectively. Since \S + M and + Q are maximally monotone on K, the operators 



B 



1 



S + M] and a := F 



-1^1 



-S + Q 



(2.17) 



are maximally monotone on /Cy. Moreover, since V is self-adjoint and p-strongly 
positive, one can easily see that weak and strong convergence in /Cy are equivalent 
with weak and strong convergence in /C, respectively. 
Consequently, for every n > we have 



V{xn - G -S + M (y„ - fo„) + -S6„ - 



1 



^ Vxn ^[V + -S + M'j (y„ - K) + [^-S + K - hi 

^XnG (id + + m)) (y„ - b„) + + 6„ - b'^ 



^y^ = (Id + B)-^ (xn - e^) + K 



(2.18) 



and 



V{Wn - Zn) G + (Z„ - dn) + ^Sdn - < 



^ Zn= [Id + V-' {-S + Q 



1 



Wn-e^) +d. 



<^ Zn = (Id + A) ^ (Wn - e^j + dn. 

Thus, the iterative rules in (2.14) become 

yn = Jb (xn - e^) + bn 

(Vn > 0) ' 



(2.19) 



Zn — J A i^Vn "^n ^nj ~^ ' 
Xn+1 = Xn + \n{Zn — Vn) 

In addition, we have 

zer {A + B)= zer (f^^ (M + 5 + Q)) = zer (M + 5 + Q) . 



(2.20) 
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By defining for every n > 

= Jb {xn - e^) -Jb {xn) + bn and cxn = J A {^Vn - Xn - s'^-Ja (2y„ - a3„) + d„, 
the iterative scheme (2.20) becomes 



(Vn > 0) 



Vn = JB{Xn) + ^n 

Zn = J A (22/„, - Xn) + OLn 

Xn-\-l — Xji + Ayi(Zji Uri) 



(2.21) 



-CXD 



thus, it has the structure of an error-tolerant Douglas-Rachford algorithm (see 
(i) The assumptions made on the error sequences yield 

+00 +00 +00 +00 

X! -^nMnlk < +00, ^ AnlKllx; < +00, ^ A„||&„||a: < +00, ^ A„||bJ^||A: < +C 
n=0 ra=0 n=0 n=0 

(2.22) 

and, by the boundedness of V^^ , S and V, if follows 

+00 +00 



^ An||e^||7C < +00 and ^ A„||e5i[||7c < +oo. 

n=0 n=0 



(2.23) 



Further, by making use of the nonexpansiveness of the resolvents, the error sequences 
satisfy 



n=0 



+00 

XI [\\0Ln\\K + IIAnIk] < ^ K \\Ja {^Vn - Xn - - Ja (2y„ - Xn) \\k + ll^nlk 

n=0 

+ \\Jb {xn - e^n) - Jb {Xn) \\k + ll^'nlk 



+00 

<EAn 
n=0 



l^nlk + Mnik + ll^nlk + ll^'nlk 



< +0O. 



By the linearity and boundedness of V it follows that 



+00 

51 [\\0Ln\\Kv + WuWkv] < +°^- 
n=0 

(i)(a) According to [9, Theorem 2.1(i)(a)] the sequence {xn)n>o converges weakly in 
JCy and, consequently, in /C to a point x S (RaRb) with Jbx G zer(A + B). The 
claim follows by identifying Jbx and by noting (2.8). 

(i)(b) According to [9, Theorem 2.1(i)(b)] it follows that {RaRbXu — Xn) — >■ (n — )■ 
+00). From (2.21) it follows that for every n > 

K{Zn - Vn) = Y (^A{RB{Xn) + 2^„) - Xn + 2an) , 
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thus, by taking into consideration the nonexpansiveness of the reflected resolvent and 
the boundedness of (A„)„>o, it yields 



\>^n{Zn - yn)\\Kv <^\\RARBXn - Xn\\Kv 



A 
2 

A 
2 

<\\RARBXn - Xn\\Kv + ^ri [ll^nlkv + ll/^nlkv] • 



+ -l^\\RA{RBXn + 2/3 J - RAiRsXn) + 2Q!„,||a:^, 



The claim follows by taking into account that A^ [lloinlkv + ll/^nlkvl — ^ (n — )■ +oo). 

(i) (c) As shown in (a), we have that a;„ — >■ S G fix {RaRb) {n — >■ +oo) with Jbx G 
zer(A+S) = zei{M+S+Q). Moreover, by the assumptions we have b„ — >■ (n — t- +oo) 
(cf. (2.13), thus ^ (n ^ +oo) (cf. (2.15)) and ^ (n ^ +oo). Hence, by the 
continuity of Jb and (2.21), we have 

= Jb {xn) + /9„ ^ Jbx G zer (M + 5 + Q) (n +oo). 

(ii) The assumptions made on the error sequences yield 

+00 +00 +00 +00 

X! Il'^^lk < +00, ^ IKIk < +00, ^ \\bn\\K < +00, WKWk < +00, 

n=0 n=0 n=0 n=0 



thus. 



This implies that 



+00 +00 

^ lle^lk < +00 and ^ \\ei\\K < +oo. 

n=0 n=0 



J2 [ll"™|k + ll^nlk] < +00 
n=0 



which, due to the linearity and boundedness of V, yields 

iW^nWKv + ll^nlkv] < 

n=0 

Since A and , i = 1, m, are uniformly monotone, there exist increasing functions 
4>A '■ — >■ [0,+oo] and : ^+ — >■ [0,+oo],i = l,...,m, vanishing only at 0, such 

i 

that 



{x-y,u- z) > (pAiWx - y\\n) V(x,n),(y,z) G gra^ 

{v-w,p-q) > (p^-i {\\v-w\\g^) \/ {v,p),{w,q) G gra^r^ Vi = l,...,m 



(2.24) 



The function (j)M '■ IR+ — ^ [0, +oo] 



9M 



(c) = inf < (pA{a) + J2'I^B7''ibi) : 



i=l 



\ 



i=l 
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is increasing and vanishes only at and it fulfills for each {x, u), (y, z) E gra M 



{x-y,u- z),^>(j)M{\\x-y\\]c). (2.26) 

Thus, M is uniformly monotone on /C. 

The function (ps '■ ^+ — ^ [0, +oo], (pB{t) = 0M (^^^7|^^^ ' ^® increasing and vanishes 

only at 0. Let be {x,u), {y, z) G gra-B. Then there exist v G Mx and w G My 
fulfining Vu = ^Sx + v and Vz = ^Sy + w and it holds 

{x-y,u- z)^^ = {x-y,Vu- Vz)^ 



x-y, Qsa; + - Qsy + w 

(2.26) 

> 4>M - y\\K) 



X - y\\K. 



K 




(2.27) 

Consequently, B is uniformly monotone on K,v and, according to [9, Theorem 
2.1(ii)(b)], {JBXn)n>o Converges strongly to the unique element y £ zer(A + B) = 
zer (M + S + Q). In the Ught of (2.21) and using that /3„ (n +oo), it follows 
that y^ ^ y (n ^ +oo). □ 

Remark 2.1. Some remarks concerning Algorithm 2.1 and Theorem 2.1 are in order. 

(i) Algorithm 2.1 is a fully decomposable iterative method, as each of the opera- 
tors occurring in Problem 2.1 is processed individually. Moreover, a considerable 
number of steps in (2.4) can be executed in parallel. 

(ii) The proof of Theorem 2.1, which states the convergence of Algorithm 2.1, relies on 
the reformulation of the iterative scheme as an inexact Douglas-Rachford method 
in a specific real Hilbert space. For the use of a similar technique in the context 
of a forward-backward- type method we refer to [14]. 

(iii) We would like to notice that the assumption X^n:^ '^'^ll'^nll'H < +oo does not 
necessarily imply neither that (||ara||-H)ra>o is summable nor that (an)n>o (weakly 
or strongly) converges to as n — ^ -|-oo. We refer to [9, Remark 2.2(iii)] for further 
considerations on the conditions imposed on the error sequences in Theorem 2.1. 

Remark 2.2. In the following we emphasize the relations between the proposed algo- 
rithm and other existent primal-dual iterative schemes. 

(i) Other iterative methods for solving the primal-dual monotone inclusion pair in- 
troduced in Problem 2.1 were given in [10] and [14] for , i = 1, m monotone 
Lipschitzian and cocoercive operators, respectively. Different to the approach 
proposed in this subsection, there, the operators D^^,i = 1, ...,m, are processed 
within some forward steps. 
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(ii) When for every i = 1, ...,m one takes Di{0) = Qi and Di{v) = G Gi \ {0}, 
the algorithms proposed in [10, Theorem 3.1] (see, also, [6, Theorem 3.1] for the 
case m = 1) and [14, Theorem 3.1] applied to Problem 2.1 differ from Algorithm 
2.1. 

(iii) When solving the particular case of a primal-dual pair of convex optimization 
problems discussed in Example 2.1 one can make use of the iterative schemes 
provided in [11, Algorithm 3.1] and [7, Algorithm 1]. Let us notice that particu- 
larizing Algorithm 2.1 to this framework gives rise to a numerical scheme different 
to the ones in the mentioned literature. 



2.2 A second primal-dual algorithm 

In Algorithm 2.1 each operator Li and its adjoint L* ,i = 1, ...,m are processed two times, 
however, for large-scale optimization problems these matrix-vector multiplications may 
be expensive compared with the computation of the resolvents of the operators A, B^^ 
and D^^ ,i = 1, m. 

The second primal-dual algorithm we propose for solving the monotone inclusions 
in Problem 2.1 has the particularity that it evaluates each operator Lj and its adjoint 
L*,i = l, m, only once. 

Algorithm 2.2. 

Let xq £ n, {yifi, . . . , ym,o) £ Gi x ■ ■ ■ x Gm, (vifi, Vm,o) & Gi X ■ ■ . X Gm, and r and 
(Tj,i = 1, ...,m, be strictly positive real numbers such that 



■^ai\\L 



III ^ 



1=1 



Furthermore, set = (j~^t X^ilLi Cill-^dP' ^ = l,...,m, let (An)n>o be a sequence in 
(0, 2), (a„)„>o a sequence in "H, (6i,n)n>o and (di,n)n>o sequences in Gi for all i = 1, . . . , m 
and set 

Pl,n = JtA {Xn - T {J2iLl L*Vi,n - z)) + ttn 

For i = 1, . . . ,m 

(Vn > 0) P2,i,n = J'y.Di {yi,n + liVi,n) + di^n 

yi,n+l = yi,n + K{P2,i,n " yi,n) 

P3,i,n = 'h B-^ (^«:« + (-^i(2pi,n " Xn) - {2p2,i,n " yi,n) " ri)) + bi^n 

i 

Vi,n+1 = Vi^n + K{P3,i,n " ^^i.n)- 

(2.28) 



Theorem 2.2. In Problem 2.1 suppose that 

/ m 

zeiB.n[A + Y,L*{BiUDi){Li 



(2.29) 



1=1 



and consider the sequences generated by Algorithm 2.2. 
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(i) If 

+00 +00 

X! -^nllanllw < +00, X! '^ri{\\di^n\\g, + \\bi,n\\g,) < +00, i = l,...,m, 

n=0 n=0 

and J2n=o An (2 - An) = +00, then 

(a) (xn, yi,n, ym,n, vi^n, • • • , fm,n)n>o Converges weakly to a point 

{x,yi,...,y^,Vi,...,Vm) G Ti X Gi X . . . X Grn X Gl X ■ ■ ■ X Gm SUch 

that {x,vi, ...,Vm) is a primal-dual solution to Problem 2.1. 

(h) \n{Pl,n-Xn) (n +Oo), Xn{P2,i,n-yi,n) (n +Oo) and Xn{p3,i,n- 

Vi,n) — )■ (n — >■ +oo) for i = 1, ...,m. 
(c) whenever % and Gi,i = l,...,m, are finite- dimensional Hilbert spaces, 
{xn, vi^n, • • • ) Vm,n)n>o convcrgcs strongly to a primal-dual solution of Prob- 
lem 2.1. 

(ii) If 

+0C +00 

y2\\an\\n < +00, 'S2{\\di^n\\gi + \\bi,n\\gi) <+oo,i = l,...,m, inf An > 

n=0 n=Q 

and A, and Di,i = 1, ...,m, are uniformly monotone, 

then {pi,mP3,i,n, ■ ■ ■ ,P3,m,n)n>o convcrgcs strongly to the unique primal-dual solu- 
tion of Problem 2.1. 

Proof. We let ^ = ^1 X ... X Q^^ be the real Hilbert space endowed with the inner 
product and associated norm defined in (2.6) and consider 

K = nxg xg, 

the real Hilbert space endowed with inner product and associated norm defined for 
X = {x,y,v),u = {u, q,p) e /C as 



{x,u)^ = {x,u)y^ + {y,q)g + {v,p)g and ||a;||K; = ^ + \\y\\g + \\v\\^, (2.30) 
respectively. In what follows we set 

y = {yi,-,ym), V = {vi,...,Vm), y = {yi,...,ym), V = {vi,...,Vm)- 

Consider the set-valued operator 
M : /C 2 , {x,y,v) {-z + Ax, Diyi, . . . , Dmym,ri + Bi + B 

which is maximally monotone, since A, Bi and Di,i = 1, ...,m, are maximally monotone 
(cf. [1, Proposition 20.22 and Proposition 20.23]) and the bounded linear operator 

S-.K^K, {x,y,v) ^ {^Y^Jj*Vi,-vi,...,-Vm-,-Lix + yi,...,-LraX + yJ^ , 
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which proves to be skew (i. e. S* = —S) and hence maximally monotone (cf. [1, Example 
20.30]). Since domS^ = /C, the sum M + S is maximally monotone, as well (cf. [1, 
Corollary 24.4(i)]). Further, we have 



i=l 



(2.29) ^ {3 X £ n) z £ Ax + ^ L* {Bi □ A) {Ux - ri) 
{3{x,v) enxQ) 
^ {3{x,v) enxQ) 



Vi £ (BiDDi) {LiX - ri) , i = 1, . . . ,m 



LiX - Vi £ B~^Vi + D~^Vi, i = 1,... ,m 



oe-z + Ax + YT=i ^v^ 

<^ (3 (x, y,v) £K) { G DiUi - Vi, i = 1, . . . ,m 

G r j + -B^^Vj - LiX + Ui, i = 1, . . . ,m 

^ (3 (x, y, v)eK) {0,...,0)e{M + S) {x, y, v) 

4^ zer (M + 5) / 0. (2.31) 



From the above calculations it follows that 
(x, y, v) G zer (M + S) ^ 



^ z-ET=iL*v^^Ax 

I Vi e {BiODi) {LiX - ri) ,i = 1, ...,m 

44> {x,vi, ...,Vm) is a primal-dual solution to Problem 2.1. 

(2.32) 



Finally, we introduce the bounded linear operator 
V :K^K 



{x,y,v) ^ i- - Y^ L*iVi, — + vi 



i=l 



71 



Um f 1 

L, . . . , \-Vm, Lix + yi, 



In 



0"! 



which is self-adjoint, i. e. V* = V. Furthermore, the operator V is p-strongly positive 
for 



1-2 



\ 



i=l 



mm 



1 1 



1 



, , . . . , , , . . . , 

T 71 7m CTl (Tn 



which is a positive real number due to the assumption 



T > CTiii^ii 
i=l 



< 



4 



(2.33) 



made in Algorithm 2.2. Indeed, for ji = SE^i ^''ill-^ilP it yields for each i 

1, . . . ,m. 



2{LiX -yi,Vi)g^ < 2||Li||||x||^||vi||g^ + 2||yj||gJ|wj|b^ 



\ i=i 



.ii2Mki 
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and, consequently, for each x = (x, v) E /C, it follows that 



+ E 



i=l L 



W^Wq. I 

7i o-j 



> 1-2, 



mm 



- 2^ (LjX - yi, 

i=l 

11 11 



The algorithmic scheme (2.28) is equivalent to 



, , . . . , , , . . . , 



T 71 7.m 0"! 



\x\ 



K 



p\\x\\k.- 



(2.34) 



(Vn > 0) 



For i = 1 , . . . , m 

+ V^,n G A(p2,i,n - di^n 
yi,n+l = yi,n + K{P2,i,n - yi,n) 

Li{Xn - Pl,n) + yi,n " P2,i,n 



(2.35) 



We introduce for every n > the following notations: 



•^n — (^ri) 2/1,71) • • • ) ym,m V\^ni • ■ ■ i Vm,n) 

,niP2,l,m ■ ■ ■ iP2,m,n-iP'i,l,m ■ ■ ■ ■,PZ,m,n) 
— (^m d\ ni • • • , dm,ni ^l,n) • • • ) &m,n) 



(2.36) 



a: 



V T 5 ^-^ , • • • 



In 



Hence, the scheme (2.35) can equivalently be written in the form 



(Vn > 0) 



V{xn - Pn) e{S + M) - an) + 5a„ - al 



(2.37) 



Considering again the Hilbert space JCy with inner product and norm respectively 
defined as in (2.16), since V is self-adjoint and p-strongly positive, weak and strong 
convergence in Ky are equivalent with weak and strong convergence in /C, respectively. 
Moreover, A = V^^ {S + M) is maximally monotone on ICy- Thus, by denoting 
e„ = {{S + V) an — a^) for every n > the iterative scheme (2.37) becomes 



(Vn > 0) 



Pn — "^A G-n) + Otn 

— "^n ~l~ ■^ni.Pn •^n) 



(2.38) 



Furthermore, introducing the maximal monotone operator B : K, ^ 2^, x i-^ {0}, 
and defining for every n > 



J A {Xn 



J A {,Xn) ~\~ an, 
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the iterative scheme (2.38) becomes (notice that Jb = Id) 

(Vn > 0) 



Vn = JB{Xn) 
Pn = J A (2y„ - Xn) + Otn 
Xn+l = Xn + ^n{Pn ~ Vn) 



(2.39) 



thus, it has the structure of the error-tolerant Douglas-Rachford algorithm from [9]. 
Obviously, zer(A + B) = zer(M + S). 

(i) The assumptions made on the error sequences yield 



X! ^nl|ctnlk < +00 and ^ A„||e„||x; < + 



oo, 



n=0 n=0 

thus, by the nonexpansiveness of the resolvent of A, 



oo 



n=0 



and, consequently, by the linearity and boundedness of V, 

+00 

X! -^n||o^nlkv < +°°- 

(i)(a) Follows directly from [9, Theorem 2.1(i)(a)] by using that Jb = Id and relation 
(2.32). 

(i)(b) Follows in analogy to the proof of Theorem 2.1(i)(b). 

(i) (c) Follows from Theorem 2.2(i)(a). 

(ii) The iterative scheme (2.38) can be also formulated as 



(Vn > 0) 

with the error sequence fulfilling 



Pn — JAiXn) + OL-n 

y„ = Jb (2p„ - Xn) , (2.40) 

Xn+l = Xn + XniVn " Pn) 



+00 

ll^nlkv < 

n=0 

The statement follows from [9, Theorem 2.1(ii)(b)] by taking into consideration the 
uniform monotonicity of A and relation (2.32). □ 

Remark 2.3. When for every i = 1, ...,m one takes Di(0) = Qi and Di{v) = Vw G 
Qi \ {0}, and [di^n)n>o as a sequence of zeros, one can show that the assertions made in 
Theorem 2.2 hold true for step length parameters satisfying 

m 

T^aiWLiW^ < 1, 
1=1 

when choosing (yi,Oi • • • , ym,o) = (0, ...,0) in Algorithm 2.2, since the sequences 
(yi ,ni • • • 5 Vm,n\i>o and {yi^ni ■ ■ ■ 1 Vm,n\i>G Vanish in this particular situation. 
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Remark 2.4. In the following we emphasize the relations between Algorithm 2.2 and 
other existent primal-dual iterative schemes. 

(i) When for every ^ = 1, ...,m one takes Di{0) = Qi and Di{v) = \/v £ and 
('^j,n)n>o as a sequence of zeros, Algorithm 2.2 with {yifl, . . . , ymfl) = (0, . . . , 0) as 
initial choice provides an iterative scheme which is identical to the one in [14, Eq. 
(3.3)], but differs from the one in [10, Theorem 3.1] (see, also, [6, Theorem 3.1] 
for the case m = 1) when the latter are applied to Problem 2.1. 

(ii) When solving the particular case of a primal-dual pair of convex optimization 
problems discussed in Example 2.1 and when considering as initial choice yi^o = 
0, Algorithm 2.2 gives rise to an iterative scheme which is equivalent to [11, 
Algorithm 3.1]. In addition, under the assumption of exact implementations, 
the method in Algorithm 2.2 equals the one in [7, Algorithm 1], our choice of 
(An)n>o to be variable in the interval (0, 2), however, relaxes the assumption in [7] 
that (A„)n>o is a constant sequence in (0, 1]. 



3 Application to convex minimization problems 

In this section we particularize the two iterative schemes introduced and investigated in 
this paper in the context of solving a primal-dual pair of convex optimization problems. 
To this end we consider the following problem. 

Problem 3.1. Let Ti and Qi, i = l,...,m, be given real Hilbert spaces, / G r('H) and 
z e 7i, gi,li G r(^j), Vi e Gi, i = l,...,m, and Li : H ^ Gi, i = l,...,m, nonzero 
bounded linear operators. Consider the convex optimization problem 

(P) mf |/(x)+£(5^n4)(L,x-r,)-(x,z)| (3.1) 

and its conjugate dual problem 

{/ m \ m 

-fiz-Y^ L*v, - J2 i9Kv^) + l*M + {v^, u)) . (3.2) 
V 1=1 / i=l J 

Considering the maximal monotone operators 

A = df, Bi = dgi and Di = dk, i = 1, m, 
the monotone inclusion problem (2.1) reads 

m 

find X such that z G df{x) + 'Y^L*{dgindli){LiX — rj, (3.3) 

i=l 

while the dual inclusion problem (2.2) reads 

find 771 G ^1 TJ^ G such that (3x en){^ ~JF=X f ^^^"^l • i 

' yvi e{dgindli){LiX - Vi), i = l,...,m. 

(3.4) 
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If (x, t"!, . . . , Vm) £ Ti X Qi . . . X Qrn IS & primal-dual solution to (3.3)-(3.4), namely, 

(3.5) 



m 



z - ^L*Vi £ df{x) and Vi G {dgiUdli){LiX - ri), i = 1, . . , 

i=l 



then X is an optimal solution to (P), (vi, . . . ,Vm) is an optimal solution to {D) and 
the optimal objective values of the two problems, which we denote by v{P) and v{D), 
respectively, coincide (thus, strong duality holds). 

Combining this statement with Algorithm 2.1 and Theorem 2.1 give rise to the 
following iterative scheme and corresponding convergence results for the primal-dual 
pair of optimization problems (P) — (D). We also use that the subdifferential of a 
uniformly convex function is uniformly monotone (cf. [1, Example 22.3(iii)]). 

Algorithm 3.1. 

Let xq £ Ti, (f 1,0, • • • , Vm,o) £ GiX . . .X and r and fij, i = 1, m, be strictly positive 
real numbers such that 

m 

T^a^\\Lif < 4. 

i=l 

Furthermore, let (A„)n>o be a sequence in (0, 2), (a„)„>o a sequence in Ti, {bi^n)n>o and 
{di^n)n>o sequences in Gi for all i = 1, . . . , m and set 

Pl^n = PrOX^j {Xn - I J2iLl L*Vi^n + Tz) + a„ 
Wl,n = 2pi,„ - Xn 

For i = 1, . . . ,m 

P2,i,n = PrOXo-.g* {Vi^n + ^LiWi^n - CTiVi) + bi^n 
W2,i,n = 2p2,i,n " ^i.n 

Xn+1 = Xn -\- An(zi^n ~ Pl,n) 



(Vn > 0) 



r v^m r*„,, 
l^i=l W2,i,', 



(3.6) 



For i = 1, . . . ,m 



Z2,i,n = PrOX^^j* {W2,i,n + ^Li{2zi^n " Wl,n)) + di^r. 
Vi,n+l = Vi^n + ^n{z2,i,n — P2,i,n)- 



Theorem 3.1. For Problem 3.1 suppose that 

z e ran [df + L*{dg,Udk){Li ■ -r,) j , 

and consider the sequences generated by Algorithm 3.1. 
(i) If 

+ 00 +00 

A„||a„||^ < +0O, K{\\di,n\\g, + ||^j,n||gj < +oo,i = 1, . . . ,m, 

n=0 n=0 



(3.7) 



and J2n=o -^n(2 - A„) = +00, then 
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(i) {xni vi^n, • • • J Vm,n)n>o Converges weakly to a point {x, vi, . . . , Vm) £ H x Qi x 
... X Qm such that, when setting 

and P2^i = Prox^^g. (^i + ^Li(2p^ -x)- air^ i = 1, m, 

Pi is an optimal solution to (P), i) • • ■ )P2 m) ^'-^ optimal solution to 
{D) and v{P) = v{D). 
(a) A„(zi,„ - pi,n) (n +00) and A„(2;2,i,„ - P2,i,n) (n +00) for 
i = I, m. 

(Hi) whenever % and Gi,i = 1, ...,m, are finite- dimensional Hilbert spaces, an — >■ 
(n — >■ +00) anc? — >■ (n — >■ +00) for i = l,...,m, then {pi^n)n>o 
converges strongly to an optimal solution to (P) and {p2,i,n, ■ ■ ■ TP2,m,n)n>o 
converges strongly to an optimal solution to (D) . 



(ii) If 

+00 +00 

llanllw < +00, V(||(ij„||g^ + ||6j„||gJ < +(X),i = 1, . . . ,m, inf A„ > 

n=0 n=0 

and f and g*,i = 1, ...,m, are uniformly convex, 

then (pi,n)n>o converges strongly to an optimal solution to (P), 
{P2,i,n, ■ ■ ■ ■,P2,m,n)n>o Converges strongly to an optimal solution to {D) and 
v{P)=v{D). 

Algorithm 2.2 and Theorem 2.2 give rise to the following iterative scheme and 
corresponding convergence results for the primal-dual pair of optimization problems 
{P)-{D). 

Algorithm 3.2. 

Let xq £ n, (2/1,0, • • • , ymfi) & Qi X . . . X Qrn, {vifl, Vmfi) & Gi X . . . X Grn, and T and 
C7i, i = 1, ...,m, be strictly positive real numbers such that 



r 



< 



1 



1=1 

-l_v^m „.|ir l|2 



Furthermore, set 7^ = TJ2iLiO'i\\Li\\ , i = l,...,m, (An)n>o he a sequence in (0,2), 
(on)n>o a sequence in H, {hi^n)n>o and ((ij_„)„,>o sequences in Qi for all i = 1, . . . , m and 
set 

pi,„ = Prox^j {xn - T {J2ill L*Vi^n " z)) + a„ 

Xn+l = Xn + Xnipi ,n Xn j 

For i = 1, . . . ,m 

(Vn > 0) P2,i,n = Prox-y-;. (yi,„ + jm^n) + di^n 

yi,n+l = yi,n + K{P2,i,n " yi,n) 

P3,i,n = PrOX<^^g. {Vi^n + CTi (Li(2pi,„ - X„) - (2p2,i,n - yi,n) " ri)) + bi^n 
Vi,n+1 = Vi^n + K{P3,i,n " Vi,n)- 

(3.8) 
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Theorem 3.2. In Problem 3.1 suppose that 

z G ran [df + ^ L* {dgi □ 9/^) (L^ • -r,) j , (3.9) 

and consider the sequences generated by Algorithm 3.2. 

(i) If 

+ 0O +CO 

5Z '^"llflnllw < +00, X! -^"(II^'.^IIg, + ll^j.nlbJ < +oo,i = 1, . . . ,m, 

n=0 n=0 

and Y.n=o An (2 - An) = +00, then 

( a) (xn, 2/i,n, • • • , ym,n, fi.n, • • • , « m,n)n>o convcrgcs Weakly to a point 
{x,yi,...,y^,vi,...,Vm) G 7^ x ^1 x . . . x x ^1 x . . . x such 
that such that x is an optimal solution to (P), (vi, ...,Vm) is an optimal 
solution to {D) and v(P) = v{D). 

(b) Xn{Pl,n-Xn) ^ {n ^ +Oo), A„(p2,i,n -^i.n) (n +Oo) and Xn{p3,i,n- 

Vi,n) — )■ (n — ^ +oo) for i = I, m. 

(c) whenever % and Qi, i = l,...,m, are finite- dimensional Hilbert 
spaces, (xn)n>o converges strongly to an optimal solution to (P) and 
(fi^n, • • • ,^^m,n)n>o Converges strongly to an optimal solution to (D). 

(ii) If 

+00 +00 

y2 \\an\\n < +00, ^{\\di^n\\gi + \\bi,n\\gi) <+oo, 1 = 1,..., m, inf An > 
n=0 n=0 

and f, li and g*,i = l,...,m, are uniformly convex, 

then, (pi,n)n>o converges strongly to the unique optimal solution to (P), 
iP3,i,n, ■ ■ ■ ,P3,m,n)n>o convcrges strongly to the unique optimal solution of (D) 
and v{P) = v{D). 

Remark 3.1. According to Remark 2.3, when U : Qi ^ M, li = (5{o}, and ((ij,n)n>o is 
chosen as a sequence of zeros for every i = 1, m, the assertions made in Theorem 3.2 
hold true for step length parameters satisfying 

m 

T^aiWUf < 1 

when taking in Algorithm 3.2 as initial choice (yi,0) ■ ■ ■ , Umfi) = (0, . . . , 0). In this case 
the sequences (?/i,„, . . . ,ym,n)n>o and (p2,i,n, • • • ,P2,m,n)n>o vanish and (3.8) reduces to 

Pl,n = PrOX^/ {Xn - T (EE^l LiVi,n - z)) + ttn 
— ~l~ An(pi,n ^n) 

Foi i = l,...,m (3.10) 

P3,i,n = PrOX^^g. {Vi^n + CFi (Lj(2pi,n - X„) - Tj)) + 6j,n 
Vi,n+1 = Vi,n + Xn{p3,i,n " Vi^n) ■ 



(Vn > 0) 
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4 Numerical experiments 



In this section we emphasize the performances of the algorithms introduced in this 
article in the context of two numerical experiments on location and image deblurring 
problems. 

4.1 The generalized Heron problem 

We start by considering the generalized Heron problem which has been recently inves- 
tigated in [12, 13] and where for its solving subgradient-type methods have been used. 

While the classical Heron problem concerns the finding of a point u on a given 
straight line in the plane such that the sum of distances from u to given points , u'^ is 
minimal, the problem that we address here aims to find a point in a closed convex set 
Q C M" which minimizes the sum of the distances to given convex closed sets 0,i C M", 
i = 1, . . . ,m. 

The distance from a point x £ M" to a nonempty set 0, C M" is given by 
d{x;n) = (II • \\n5n){x) = inf ||x-z||. 

Thus the generalized Heron problem reads 

m 

inf Vd(x;17,), (4.1) 
x^Q — f 

1=1 

where the sets J7 C M" and C M", i = 1, . . . ,m, are nonempty, closed and convex. 
We observe that (4.1) perfectly fits into the framework considered in Problem 3.1 when 
setting 

/ = 5n, and gij = || • ||, 1^ = Sq. for alH = 1, . . . ,m. (4.2) 

However, note that (4.1) cannot be solved via the primal-dual methods in [10] and [14] 
since they require the presence of at least one strongly convex function (cf. Baillon- 
Haddad Theorem, [1, Corollary 18.16]) in each of the infimal convolutions || • ||n(5Q., 
i = 1, . . . ,m, fact which is obviously not the case. Notice that 

g* : M, g*{p) = sup {{p,x) - \\x\\} = 5b(o,i)(p), ^ = 1, ■■■,m, 

thus the proximal points of /, g* and l^,i = 1, m, can be calculated via projections, 
in case of the latter via Moreau's decomposition formula. 

We test our algorithms on some examples taken from [12, 13] and denote in the final 
evaluations for every A; > by Vfc the value of the objective function of (4.1) at each 
iterate. 

Example 4.1 (Example 5.5 in [13]). Consider problem (4.1) with the constraint set Q 
being the closed ball centered at (5, 0) having radius 2 and the sets ilj, i = 1, . . . , 8, 
being pairwise disjoint squares in right position in (i. e. the edges are parallel to 
the X- and y-axes, respectively), with centers (—2,4), (—1,-8), (0,0), (0,6), (5,-6), 
(8, —8), (8, 9) and (9, —5) and side length 1, respectively (see Figure 4.1). 
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When solving this problem with Algorithm 3.1 and Algorithm 3.2 and the choices 
made in (4.2), the following formulae for the proximal points involved in their formula- 
tions are necessary: 

Froxrf (x) = (5, 0) + argmin J \\y - {x - (5, 0)) f = (5,0) + ^'5(0,2) (x - (5, 0)) 

yeB{0,2) ^ 

Pi'oxCT^g.(p) = argmin ^llz - = Pb(o,i) {v) 
z£B{0,l) ^ 



,(1.3) 



Prox,^^. (p) = p- ai Prox^-i^^ 



p — ai arg mm - 



Hence, when choosing xq £ and (vi,0; • • • , ^^8,o) G 



= p- (TiVn, ( — 

C,; \(Ti 

X ... X as starting values, 

letting (A„)„>o ^ (0,2) and r, fij G i = 1,...,8, such that TX]f=ifj < 4, the 

iterative scheme Algorithm 3.1 becomes 

Pl,n = (5, 0) + Vb{Q,2) {xn - § Ei=l Ui,n - (5, 

For i = 1,...,8 

P2,i,n = ^B(0,1) (^i,n + f ^1.") 
■l«2,i,n = 2p2 



(Vn > 0) 



Xn+l = Xn + Xn{zi^n " Pl,n) 

For i = 1,...,8 

^'j.n+l = Vi^n + An(^2,j,n — P2,i,n)- 



For the same initial choices, but when rj^f^i fij < ^, the iterative scheme in Algorithm 
3.2 reads 

Pl,n = (5, 0) + Pb(o,2) - T X;f=i - (5,0)) 

Xn+l — Xn ~\~ ^n(,Pl,n Xn") 

For i = 1,...,8 

(Vn > 0) p2,i,n = T'q, (yi,n + liVi,n) 

yi,n+l = yi,n + K{P2,i,n - yi,n) 

P3,i,n = 'Pb(0,1) {Vi,n + CTi ((2pi,rt - Xn) - (2p2,j,n " yi,n))) 
Vi,n+1 = Vi,n + K{P3,i,n " Vi,n) ■ 



k 


Algorithm 3.1 


Algorithm 3.2 


Pi,k 


Vk 


Pi,k 


Vk 





(5,-2) 


54.418914 


(5,-2) 


54.418914 


5 


(3.344027,-1.121496) 


53.046330 


(3.809999,-1.607451) 


53.174978 


10 


(3.389398,-1.185733) 


53.043638 


(3.441673,-1.253641) 


53.046054 


20 


(3.392361,-1.189747) 


53.043627 


(3.392712,-1.190221) 


53.043627 


50 


(3.392688,-1.190188) 


53.043627 


(3.392688,-1.190188) 


53.043627 


10^ 


(3.392688,-1.190188) 


53.043627 


(3.392688,-1.190188) 


53.043627 
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10r 
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- 0- 
-2 
-4- 
-6 
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Initializations 


Algorithm 3.1 


ai = 0.5, T = 0.24, 
A„ = 1.8, .To-(5,2), 
Vifl = 0, i 1, . . . ,8 


Algorithm 3.2 


Oi = 0.1, T = 0.24, 
A„ = 1.8, xo = (5,2), 
Vi.o = 0, i = 1, . . . ,8 



Figure 4.1: Example 4.1. Generalized Heron problem with squares and disc constraint set on 
the left hand side, initializations for the iterative schemes on the right hand side. 



As pointed out in the above table, our methods need less than 50 iterations to obtain 
an accuracy up to 6 decimal places for both function values and iterates. By comparison, 
the subgradient-type method in [13], when applied to this particular problem with the 
same starting point, requires more than a million iterations. 

Example 4.2 (Example 4.3 in [12]). In this example we solve the generalized Heron 
problem (4.1) in M'^, where the constraint set $7 is the closed ball centered at (0,2,0) 
with radius 1 and = 1,...,5, are cubes in right position with center at (0,-4,0), 
(-4,2,-3), (-3,-4,2), (-5,4,4) and (-1,8,1) and side length 2, respectively. The 
plotted and the initializations for Algorithm 3.1 and Algorithm 3.2 are shown in Figure 
4.2. 




Initializations 


Algorithm 3.1 


ai = 0.4, T = 0.99, 
A„ = 1.8, xo = (0,2,0) 
Vifl = 0, i = 1, . . . ,8 


Algorithm 3.2 


(T, = 0.05, T = 0.59, 
A„ = 1.8, = (0,2,0), 
Wi,o = 0, i 1, . . . , 8 



Figure 4.2: Example 4.2. Generalized Heron problem with cubes and ball constraint set on the 
left hand side, initializations for the applied methods on the right hand side. 
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Algorithm 3.1 


Algorithm 3.2 


k 


Pi,k 


Vk 


Pi,k 


Vk 





(0,2,0) 


24.18180 




(0,2,0) 


24.18180 


5 


(-0.92380,1.62587,0.08140) 


22.23482 


(- 


-0.93595,1.66118,0.09588) 


22.23627 


10 


(-0.92525,1.62890,0.07875) 


22.23480 


(- 


-0.92561,1.62957,0.07762) 


22.23480 


20 


(-0.92531,1.62907,0.07883) 


22.23480 


(- 


-0.92520,1.62880,0.07882) 


22.23480 


50 


(-0.92531,1.62907,0.07883) 


22.23480 


(- 


-0.92531,1.62907,0.07883) 


22.23480 


10^ 


(-0.92531,1.62907,0.07883) 


22.23480 


(- 


-0.92531,1.62907,0.07883) 


22.23480 



According to the above table, both algorithms solve this problem in few steps. Especially 
Algorithm 3.1 achieves accuracies up to 5 decimal places for both objective function 
values and iterates in less than 20 iterations. 



Example 4.3 (Example 5.4 in [13]). Consider again the generalized Heron problem 
(4.1) for squares in right position in R^, but this time subject to a real line. More 
concretely, we take Q. = {(xi, X2) € : (xi, X2) = (1, 6) -|- 0), t G M} and for i = 
1,...,5, the squares with center (—6,-9), (—5,4), (0,-7), (1,0) and (8,8) and side 
length 2. Figure 4.3 shows the plotted result and the initializations for Algorithm 3.1 
and Algorithm 3.2. 

10r 



Initializations 


Algorithm 3.1 


Gi = 0.1, T = 3.99, 
A„ = 1.7, .To = (-1,6), 
Vifi = 0, « = 1, . . . , 8 


Algorithm 3.2 


CTi = 0.1, T = 0.49, 
A„ = 1.7, xo = (-1,6), 
Vi.Q = 0, t = 1, . . . , 8 



Figure 4.3: Example 4.3. Generalized Heron problem with squares and line constraint on the 
left hand side, initializations for the applied methods on the right hand side. 



k 


Algorithm 3.1 


Algorithm 3.2 


Pi,k 


Vk 


Pi,k 


Vk 





(-1,6) 


42.883775 


(-1,6) 


42.883775 


5 


(-1.215422,6) 


42.884811 


(-1.136966,6) 


42.883775 


10 


(-1.093321,6) 


42.882115 


(-1.107478,6) 


42.882444 


20 


(-1.094633,6) 


42.882115 


(-1.094886,6) 


42.882145 


50 


(-1.094773,6) 


42.882115 


(-1.094773,6) 


42.882115 


10*^ 


(-1.094773,6) 


42.882115 


(-1.094773,6) 


42.882115 



Algorithm 3.1 and Algorithm 3.2 achieve accuracies for both the values of the objective 
function and the iterates in less than 50 iterations. By comparison, as shown in [13, 
Example 5.4], the subgradient-type methods are not able to provide these accuracies 
not even after 10 million iterations. 
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4.2 Image deblurring 

The second numerical experiment concerns the solving of an ill-conditioned linear inverse 
problem arising in image deblurring. To this end, we consider images of size M x N 
as vectors x G for n = MN, where each pixel denoted by Xij, 1 < i < M, 1 < 
j < N, ranges in the closed interval from (pure black) to 1 (pure white). For a given 
matrix A G M"^" describing a blur (or averaging) operator and a given vector b E M" 
representing the blurred and noisy image, our aim is to estimate the unknown original 
image x G M" fulfilling 

Ax = b. 

To this aim we are solving the following regularized convex nondifferentiable problem 
jnf^ {px - b\\^ + 02 \\Wx\\^ + aiTV{x) + 5[o,i]"(x)}, (4.3) 

where the regularization is done by a combination of two functionals with different 
properties. Here, ai, 02 G 1^++ are regularization parameters, TV : M" — )■ M is the 
discrete isotropic total variation function and W : M"" — )• M" is the discrete Haar wavelet 
transform with four levels transforming the image into wavelet coefficients with respect 
to the orthonormal Haar wavelet basis. Notice that the norm of the operator W is 
\\W\\ = and that none of the functions occurring in the objective function of (4.3) 
is differentiable. 

The picture undergoes a Gaussian blur of size 9x9 with standard deviation 4, as 
done in [5, Section 4.2], yielding a blurring operator A with = 1 and A* = A. For 
the discrete isotropic total variation functional 

A/-1 7V-1 



i=i j=i 

M-l N-1 

+ X! — Xi^nI + X! \^M,j+l — XM,j\ , 

i=l j=l 

where reflexive boundary conditions are assumed, it holds that TV{x) = ||Lx||x, where 
y = W X M", the operator L : M" ^ 3^, Xij ^ {LiXij, LsXjj), 



Xi+ij - Xij, ifi < M _ j Xij+i - Xij, a j < N 

0, if i = M and 7.2:Ei,, - < if j = iV ' 



represents a discretization of the gradient in horizontal and vertical direction with ||L|| < 



8 and || • ||x : M, \\{p, q)\\x = J2i=i J2j=i yPij + is a norm on the Hilbert space 

y. 

Consequently, the optimization problem (4.3) can be equivalently written as 

inf {fix) + 5i(Ax) + g2{Wx) + gsiLx)}, (4.4) 

where / : ^ 1, f{x) = 5[o,i]"(x), g, : W ^ R, g^iy) = \\y - b\\i, 52 : ^ M, 
g2{y) = a2\\y\\i and g-^ : y ^ R, g^iy^z) = oi ||(y, z) || x . For every p E we have 
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gtip) = +P^b and g^{p) = (5[_«2,a2]"(p) (see, for instance [2]), while, for every 

{p,q) ey it holds g3{p,q) = Ss{p,q), where (cf. [5]) 



S=\{p,q)€y: ^max^ ,JpIj + qfj < oi 



It is easy to see that for all x,p,q £ M" it holds 



1 



ProxT-/(x) = argmin -\\z — x|p = "Prg 
^£[0,1]" 2 



Prox<^^„*(p) = avgrnin \ aiz^b + l-\\z - p\\^\ = V^i ^^r^{p - aib) 
Fvox^^g*{p) = argmin -pf = 7^[-a2,a2]"(p) 



Z&[— 02,021 

Prox<^3(,*(p,g) = argmin ;^||(2:i,Z2) - {p,q)f = Vs{p,q) 
izi,z2)es^ 

where the projection operator Vs : 3^ — 5* is defined via 



{Pi,j,Qi,j) ^ "1- 



1 < i < M, 1 <j <N. 



max 



(Vn > 0) 



Hence, when choosing xq £ M" and (vi,o, V2,o, 1^3,0) £ x x 3^ as starting values, and 
letting (A„)n>o ^ (0, 2) and r, cJi, fT2, £ M++ be such that r (cji + (T22~-^^ + 8173) < 4, 
the iterative scheme in Algorithm 3.1 becomes 

Pl,n = 'P[0,1]" i^n - \ (^*Wl,n + W*V2^n. + -^*f3,n)) 
P2,l,n = ^[-1,1]" [Vx^n + ^^(2pi,„ - X„) - Oxh) 
P2,2,n = 'P[-a2,a2r (^2,n + ^H^(2pi,„ - X„)) 
P2,3,n = Vs {V3,n + f L*(2pi,„ - 

Zi,„ = 2pi,„ - X„ - I (A*(2p2,l,n - fl,n) + W^*(2p2,2,n " 'V2,n) 

+L*(2p2,3,n -t'S.n)) 
Xn+1 = Xn + Xn{zi,n - Pl,n) 

Vl,n+1 = Vl,n + A„ {P2,l,n " Vl,n + ^^(22;i,„ - 2pi,„ + 
l'2,n+l = V2,n + K {P2,2,n " f2,n + ^T^(2zi,„ - 2pi^n + X„)) 
t;3,„+l = f3,„ + A„ (P2,3,n - V3^n + ^L{2zi^n " 2pi,„ + Xn)) ■ 

Similarly, taking also into account Remark 3.1, Algorithm 3.2 can be also im- 
plemented to this problem, this time by choosing r, ui, cr2, 1T3 G such that 
r (CTI + £722-^^ + 8(73) < 1. 

Figure 4.4 shows the performance of Algorithm 3.1 and Algorithm 3.2 when solving 
(4.4) for «! = 3e-3, a2 = 2e-5, starting points xq = b and (t'i,o, "^2,0, "^^3,0) = (0,0,0) and 
parameters 



• DRl (Algorithm 3.1): cti = 1, (73 = 1, 0-3 = 0.05, r = 4 (tri + (J22'^^ + Scrg) - 0.01, 
A„ ~ 1.5 for every n > 0; 
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• DR2 (Algorithm 3.2): cti = 1, ctz = 0.05, erg 0.05, r = 1 (cti + a22-^^ + 80-3) - 0.01, 
A„ = 1.6 for every n > 0, 



and of the iterative scheme designed in [10, Theorem 3.1] for 



FBF ( [10, Theorem 3.1]): e 



20(^1+2^^16+8+1)' " ^1+2-16+8 

within the first 200 iterations, when applied to the 256 x 256 cameraman test image. 
Figure 4.5 shows the original, observed and reconstructed versions of the 256 x 256 
cameraman test image. 



for every n > 0, 



Function values 



ISNR values 



400 
350 
300 
250 
200 
150 
100 
50 




DR1 

■ - DR2 
- - FBF 



50 



100 
Iterations k 



150 



200 




200 



Figure 4.4: The evolution of the values of the objective function and of the ISNR (improve- 
ment in signal-to-noise ratio) for Algorithm 3.1 (DRl), Algorithm 3.2 (DR2) and the forward- 
backward-forward method (FBF) from [10, Theorem 3.1]. 



(a) Original image 



(b) Blurred and noisy image (c) Reconstructed image 




Figure 4.5: Figure (a) shows the clean 256 x 256 cameraman test image, (b) shows the image 
obtained after multiplying it with a blur operator and adding white Gaussian noise with standard 
deviation 10^'^ and (c) shows the reconstructed image generated by Algorithm 3.1. 
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